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Abstract 

Background: In Saccharomyces cerevisiae, structural bistability generates a bimodal expression of the galactose 
uptake genes (GAL) when exposed to low and high glucose concentrations. This indicates that yeast cells can decide 
between using either the limited amount of glucose or growing on galactose under changing environmental 
conditions. A crucial requirement for any plausible mechanistic model of this system is that it reproduces the 
robustness of the bistable response observed in vivo against inter-individual parametric variability and fluctuating 
environmental conditions. 

Results: We show how a control-theoretic analysis of the robustness of a model of the GAL regulatory network may 
be used to establish the model's plausibility in characterizing the persistent memory of different carbon sources, 
without the need for extensive simulations. Chemical Reaction Network Theory is used to establish that the proposed 
network model is compatible with structural bistability. The robustness of each of the two operative conditions 
against fluctuations of the species concentrations is demonstrated by studying the Domains of Attraction of the 
corresponding equilibrium points. Finally, we use a global robustness analysis method based on Semi-Definite 
Programming to evaluate the modification of the bistable steady states induced by multiple parametric variations 
throughout bounded regions of the parameter space. 

Conclusions: Our analysis provides convincing evidence for the robustness, and hence plausibility, of the GAL 
regulatory network model. The proposed workflow also demonstrates the power of analytical methods from control 
theory to provide a direct quantitative characterization of the dynamics of multistable biomolecular regulatory 
systems without recourse to extensive computer simulations. 

Keywords: Galactose network, Bistability, Robustness, Domain of attraction. Bifurcation, Local sensitivity. Global 
sensitivity 



Background 

Although yeasts, in common with most cellular organ- 
isms, can derive energy from a variety of different 
molecules, glucose is well-known to be their preferred 
source, because it provides more energy than any other 
saccharide. Therefore, yeasts have evolved a complex 
genetic network to make sure they can consume as 
much glucose as possible when it is available [1]. In [2], 
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the authors experimentally investigated the regulation of 
galactose metabolism in S. Cerevisice, which is mediated 
by several positive and negative feedback loops acting at 
the transcriptional level. To probe the system for mul- 
tistability, two identical cell populations were grown on 
different media, with and without galactose, respectively. 
In engineering terms, this amounts to initializing the sys- 
tem at two different operating conditions. Starting from 
these conditions, the two populations were then exposed 
to identical galactose concentrations for a period long 
enough to guarantee the attainment of steady-state con- 
ditions. For intermediate levels of the input (galactose 
concentrations), the two populations settled on different 
steady states, thus confirming the multistable nature of 
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the system. These and other experimental results have 
revealed that the GAL system exhibits bistable dynamics 
and that such bistability generates a persistent memory 
of the type of carbon source consumed by the cell in the 
past. 

In previous work, classical mathematical tools, such 

as bifurcation analysis, have been used to examine the 
dynamics of the GAL regulatory network, see e.g. [3,4]. 
Recently, we showed, using a control-theoretic analy- 
sis, that the GAL network simplified mass-action mod- 
els proposed in [5,6], do not reproduce the bistable 
behavior exhibited by the experimental studies of Acar 
and coworkers; this finding motivated us to propose 
a new model of the GAL system, [7]. In this paper 
we extend this model and provide a thorough charac- 
terization of its dynamical properties, with the aim of 
validating it as a plausible mechanistic explanation of 
the persistent memory property. Our approach starts 
with the analysis of the model's bistable dynamics as 
a structural property, arising from the topology of the 
reaction network. Afterwards, we focus on the study 
of the robustness of bistability both against fluctua- 
tions of the concentrations of the molecular species, 
caused by endogenous stochastic noise or by exogenous 
perturbations, and in the face of parametric uncertain- 
ties. The principle underpinning these analyses is that 
the quality of a model cannot be solely evaluated by 
its capability to reproduce a particular set of experi- 
mental measurements. Indeed, a common problem in 
modeling biological networks is that alternative, struc- 
turally different models can fit experimental data equally 
well [8]. In order to represent a plausible mechanis- 
tic description of a biological phenomenon, the model 
must also replicate an essential feature of biological sys- 
tems, that is robustness against inter-individual para- 
metric variability and in vivo fluctuating environmental 
conditions [9-12]. 

Our characterization of the robustness properties of 
the model starts with an analysis of the Domains of 
Attractions (DAs) of the bistable system. Roughly speak- 
ing, the DA of an equilibrium point Xe is a region V 
in the state space, such that Xg e V and every state 
trajectory crossing T> converges asymptotically to Xg. 
DA analysis is crucial for establishing whether the pro- 
posed model provides a plausible explanation of the phe- 
nomenon under investigation, since the system is actually 
able to operate around a given equilibrium point with 
some degree of robustness in the face of both intrin- 
sic stochastic noise and exogenous perturbations only if 
that equilibrium point possesses a nontrivial DA. Note 
that the estimation of the DA is, in general, a diffi- 
cult problem for systems of nontrivial dimension. In our 
approach we show how, for any mass-action model, it is 
possible to apply a convex optimization-based method. 



devised in a purely theoretical context by our group in 
[13,14], to test whether an assigned polytopic subset of 
the state space belongs to the DA of an equilibrium 
point. 

We next consider the robustness of the model's bistable 
dynamics in the fece of uncertain parameter values. Many 

examples can be found in the literature of studies apply- 
ing local sensitivity and bifurcations analysis as tools for 
characterizing the parametric robustness of biological sys- 
tems, e.g. [15-17]; however, these tools suffer from a 
significant limitation due to their inability to take into 
account more than one or two simultaneous parameters 
variation at the same time. Multi-parametric sensitivity 
analysis of biomodels is typically performed by resorting 
to extensive sampling of the admissible parameter space, 
[18,19], which requires a large large computational effort 
and can only provide probabilistic conclusions. To over- 
come these limitations, besides applying local sensitivity 
and bifurcations analysis, we employ a global sensitiv- 
ity analysis method proposed in [20,21]. This method 
is aimed at computing an outer approximation of the 
region of the state space that contains all the equilib- 
rium points of a given biosystem for all admissible values 
of the parameters. In our analysis, we devise a straight- 
forward way to adapt this method to provide robust- 
ness certificates for bistability in the face of parametric 
uncertainty. 

Thus, beyond our primary goal of validating a 
new model of the bistable GAL regulatory network, 
we also present what should be a widely applicable 
and effective procedure for investigating the plausi- 
bility of dynamical models of multistable biomolecu- 
lar circuits, without recourse to large-scale numerical 
simulations. 

Results 

A new model of the OAL regulatory network In 5. Cerevisise 

The regulatory network of galactose metabolism, depicted 
in Figure 1, is governed by the following factors: a tran- 
scriptional activator protein Gal4p, a signal transducer 
protein Gal3p and an inhibitor protein Gal80p. In the 
presence of galactose, Gal4p activates transcription of 
GAL2, GALS, GAL80, which are regulatory genes, and 
of GALl and several other genes (not shown in the 
figure), which encode the enzymes of the Leloir pathway 
(the GAL genes) of galactose metabolism. The protein 
encoded by gene GAL2 acts as a mediator of galactose 
transport into the yeast cell. In the absence of exter- 
nal galactose, GalSOp binds to the activation domain 
of Gal4p, thus inhibiting the expression of the GAL 
genes. In the presence of galactose, however, the inducer 
GalBp is activated to form the complex Gal80p:Gal3p', 
which promotes the shuttling of GalSOp from the nucleus 
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to the cytoplasm. This decreases the fraction of GalSOp- 
bound Gal4p in the nucleus. Thus, galactose relieves the 
inactivation of Gal4p and promotes transcription of the 
GAL genes [1]. 

The mathematical model of this regulatory network 
considered here is based on mass-action kinetics and rep- 
resents an extended version of the model proposed in 
[7]. In the new version of the model, the reaction of 
reversible dissociation of the complex Gal80p:Gal3p'' is 
explicitly included, since this was found to be essential 
to ensure robustness of the bistable dynamics, according 
to the analysis procedure that will be illustrated in the 
next sections. Moreover, in the new model also the Gallp 
protein dynamics are taken into account, since the con- 
centration of this protein is taken as a measure of Gal4p 
activity in the experiments reported in [2]. model con- 
sists of the following set of nonlinear ordinary differential 
equations (ODE's), 



G3 = /C8G4 - /c2G3Gint + kr2G3a - M13G3 (la) 
Gint = kiGexGl — /(iG^Gint + KlG'ia — MieGint (lb) 
Gsa = A'2G3Gint — kylGza — k^Gi^^G^a — f^sGsa 

+ /Cr4G80,3flG4 — krigGgoGsa + /figGgo.Sa (Ic) 

G4 = /C5 — /C11G4G8O + /CrllG4_80 + /C4G4_8oG3a 

— /Cr4G80,3flG4 — fleG^ (Id) 
Gso = — /fl 1 G4 G80 + /Crl 1 G4,80 +/<^7 G4 — /c^ig GgoG^a 

+ ki9Gso,3a~ t^iiGso (le) 
G4,80 = /<'llG4G80 — /<'rllG4,80 — /<'4G4,8oG3a 

+ /Cr4G80,3flG4 — Ml2G4,80 (If) 
G80,3fl = ^4G4,8oG3a — A'r4G80,3flG4 — /Xi5G80,3a 

— kigGsOM + krigGsoG^a (Ig) 

G2 = /C9G4 - M17G2 (Ih) 
Gi = /C10G4 - /xigGi , (li) 



Galactose 



Cytoplasm 




Nucleus 



Figure 1 Schematic diagram of the CAL regulatory network In S. Cerevisise. Galactose import is guaranteed by Gal2p; internalized galactose 
activates Gal3p, which sequesters GalSOp in the cytoplasm, shuttling it from the nucleus. The transcriptional activator Gal4p, which is constitutively 
bound to promoters of the GAL genes, is then released from the inhibitory action of GalSOp and activates expression of the GALl, GAL2, GALS and 
G/\1.80 genes. The regulatory network features two positive and one negative feedback loops. 
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where the description of each state variable is reported 
in Table 1. The total concentration of external galactose 
Gex = constant. 
The ODE model (1) can be rewritten in compact form as 

X = Nvix,p) , (2) 

where the species concentrations, namely the state 
variables listed in Table 1, are denoted by x := 
(G3GintG3aG4G8oG4,8oG8o,3flG2Gi)^, the parameters 
by p:={kik2---kri9)'^ (the full list of param- 

eters is reported in Table 2), N e M^^^^ is the sto- 
ichiometric matrix and v(x,p) e R^^ is the vector 
of reaction rates. Since the parameters are inher- 
ently positive, positive values of x will result in 
positive values of v{x,p), i.e. x € =^ v{x,p) € K^^. 
In the next section, we present the results of our 
analysis of the bistable dynamical properties of this 
model. 

Structural analysis of the proposed model's network 
topology confirms bistable dynamics 

The first step of our procedure consists in the analysis 
of the topology of the regulatory network, to determine 
whether its structure can admit a bistable behavior Sub- 
sequently, we will determine a possible realization (i.e. 
a set of parameter values) of model (1) that exhibits 
bistability. 

The persistence of cellular memory exhibited by the 
galactose regulatory network is a system-level property 
which results from the interactions of several species 
in multiple nested feedback loops. Two coupled pos- 
itive feedback loops involve the galactose permease 
Gal2p and the signaling protein Gal3p, while a nega- 
tive feedback loop involves the inhibitor GalSOp. Recall 



that, according to [22], the existence of a positive- 
feedback loop, or a mutually inhibitory, double-negative- 
feedback loop, is a necessary condition for the occurrence 
of multistability. 

Model (2) is said to exhibit bistability if there exist a 
parameter vector p e E^^, and two finite distinct equilib- 
rium points Xei,Xe2 £ such that 

Nv(Xei,p) = 0 (3a) 
Nv{xe2,p) = 0 (3b) 

The existence of a solution to Eqs. (3) can be determined 
by using Chemical Reaction Network Theory (CRNT) 
[23,24], which provides a straightforward way to establish 
whether an assigned network structure can admit multiple 
steady states. Furthermore, CRNT provides an algorithm 
to compute a feasible value of the parameter vector p and 
the associated equilibrium points Xei, and Xe2- This anal- 
ysis confirms that system (1) admits two asymptotically 
stable equilibrium points, (reported in Table 1), when the 
kinetic parameters take the values in Table 2. Actually, the 
values reported in Table 2 have been obtained by scaling 
to appropriate dimensions the numerical values returned 
by the CRNT algorithm. The scaling is required to obtain 
a good agreement with previously published experimen- 
tal data and numerical simulations, [2,25,26]. First of all, 
CRNT returns pure numbers, without associated units; in 
the light of previously published results, we have scaled 
all the quantities by a factor 10~^, such that the concen- 
trations are in the order of /iM. Subsequently, since in the 
published experiments the external galactose concentra- 
tion, Gex, lies in the order of mM, while the CRNT algo- 
rithm returned a unit value at equilibrium, we have scale 
by a factor 10~^ the kinetic constant ki, which appears in 
the galactose import term k\GexG2 in equation (lb). Such 
scaling yields an equilibrium concentration of ImM for 



Table 1 Steady states of the mass-action model (1), with the parameters values given in Table 2 





Species 


Description 


Xe2 


172.8212 


G3 


GaBp protein 


2711.1839 


172.8208 


Gint 


internalized galactose 


2711.2003 


1.0 


63<J 


active GaBp protein 


318.5443 


1.0 


G4 


Gal4p protein 


19.9479 


1.0 


Gao 


GalSOp protein 


0.3061 


21.0604 


G4,80 


Gal4p:Gal80p complex 


2.1126 


7,5945 


G803a 


Gal80p:Gal3p active complex 


589.1342 


1.0 


G2 


GaBp protein 


19.9479 


1.0 


Gi 


Gall p protein 


19.9479 


1000.0 


Gex 


external galactose 


1000.0 



Concentration values are all given as [nM\ units. 
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Table 2 A set of parameters values that renders system (1 ) bistable 





Value 


Parametsr 


Wall la 

Value 




U.lo 14 L/^/M-hJ 


/, 


or" oior" n /U1 

85.8185 [1/hJ 


h 


8.458Db-4 [l/ftM-nJ 


km 


4.7482b-2 [1//x/M-hJ 




16.6691 Ll/nJ 


Atl2 


1.0 [1/h] 




1 .0 [1/h] 


Atl3 


1.0 [1/h] 


k4 


3.0749 [V/iM-h] 


M14 


1.0 [1/h] 




0.1317 [1/h] 


/il5 


1.0 [1/h] 


ks 


22.0604 [/iM/h] 


Atl6 


1.0 [1/h] 




1.0 [1/h] 


Atl7 


1.0 [1/h] 


kj 


29.6549 [1/h] 


Mis 


1.0 [1/h] 


ks 


181.4157 [1/h] 


/t,9 


36.6342 [1//i/M-h] 


kg 


1.0 [1/h] 




222.0536 [1/h] 


km 


1.0 [1/h] 







The values have been computed through the CRNT algorithm and then scaled to suitable dimensions (see the results Section 'Structural analysis of the proposed 
model's network topology confirms bistable dynamics'). 



Gex', note also that it does not affect other equations, since 
the kinetic constant k\ does not appear elsewhere in the 
model. 

Characterization of the domains of attraction confirms 
robustness of the bistable equilibria 

Subsequently to the determination of the asymptotically 
stable equilibrium points, a primary goal in the char- 
acterization of the behavior of a system is that of 
estimating the DA's of such points. Accurate esti- 
mates of the DA's provide valuable information about 
the ability of a system to reject perturbations driv- 
ing the system away from its steady state condi- 
tion. At the same time, the boundaries of the DA's 
constitute the concentration thresholds for the activation 
of the switching mechanism between different operative 
conditions. 

The methodology proposed in [13], which allows to 
check whether an assigned box in the state space belongs 
to the DA of an equilibrium, has been employed in our 
study. It is worth noticing that the main result of [13] leads 
to a Linear JVIatrix Inequality (LMI) feasibility problem, 
which can be solved efficiently via off-the-shelf numerical 
algorithms. 

In order to find the largest possible estimates of 
the DA's of Xei and Xei, namely 2?i and "Di, our 
procedure takes two small initial polytopic regions, 
surrounding the equilibrium points, and then iteratively 
stretches them along the different dimensions of the 
state space until the feasibility conditions are no longer 
verified, thus obtaining two inner approximations of 
the DA's. 



The estimates obtained by means of this procedure are 

Vi = [ 100.82, 312.82] x [ 100.82, 332.82] x [ 0.0, 4.0] 
X [ 0.0, 4.0] X [ 0.7, 3.0] X [ 18.06, 24.06] 
X [ 2.59, 27.59] x [ 0.0, 3.0] x [ 0.0, 3.0] 

©2 = [ 211.18, 5213.2] X [ 209.20, 5622.0] 
x[ 18.54, 675.5] x[ 8.95, 30.9] 
x[ 0.0, 0.7] x[ 1.11, 3.2] 
X [ 89.13, 1528.1] X [ 6.9479, 58.9] 
x[7.95, 320.9] 

for Xe\ and Xe2, respectively (the two boxes are plotted in 
normalized parallel coordinates in Figure 2a). The validity 
of these estimates is confirmed by numerical simulations, 
performed with initial conditions varying within the boxes 
computed by the proposed approach (see Figure 3). Note 
that the admissible excursion intervals, as determined by 
the estimated DA's (reported in Figure 2a), are fairly large 
for most of the state variables: looking at Figure 2a one 
can readily recognize that the key species that drives the 
switching between the two metabolic conditions is the 
complex Gal4p:Gal80p, which is associated to a smaller 
admissible fluctuation interval with respect to the other 
species (especially in the low galactose concentration con- 
dition). Thus, the DA's analysis highlights that a tight regu- 
lation of the concentration of Gal4p:Gal80p is paramount 
to the proper functioning of the genetic switch. 
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Figure 2 Estimated regions in the state space, a) Estimated DA's of tlie two equilibrium points; b) Initial guess for the computation of the robust 
steady state subsets; c)-h) Robust steady state boxes corresponding to several admissible range of simultaneous variation of the parameters 
kt , k2, ks, kj , kg, kg, fit3, fi-\e, fiu ■ In all panels both the low (dark gray) and high (light gray) galactose concentration conditions are considered. 
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Figure 3 Free evolutions, for different initial conditions, of the concentrations of Gal4p and Gall p proteins with the set of parameters 
given by application of CRNT toolbox. The curves funnel into low (in dark gray) or high (in light gray), depending on initial conditions, confirming 
the bistable nature of the system. 



Local and global analysis confirm robustness of the 
bistable equilibria to parametric uncertainty 

In this section we provide further support for the plau- 
sibiUty of the proposed model of the GAL regulatory 
system by characterizing its robustness with respect to 
parametric uncertainties. 

The underlying principle is that, in view of the large 
inter-individual variability of biochemical parameters, for 
a model to be considered plausible it is not sufficient 
to reproduce the qualitative behavior of the biological 
system for just one set of parameter values; instead, this 
behavior must be exhibited over a nontrivial subset of the 
parameters space. 

First, a classical sensitivity analysis is performed by 
employing the method discussed in [27]: the state vari- 
ables ode's are coupled to the equations of sensitivity. 

This allows us to compute a numerical solution to 
the whole set of equations, thus simultaneously obtain- 
ing both the state variables and the associated sensitivity 
coefficients. 

The normalized sensitivity coefficients for the proposed 
model are shown in Figure 4: greater sensitivity is exhib- 
ited by the parameters involved in the feedback terms 
{ky, ks, kg), the basal expression of Gal4p (/cs), those 
involved in the internalization of external galactose and in 
the activation of Gal3p (/ci, /C2), and the parameters that 
describe the degradation of Gal3p, internalized galactose 
and Gal2p (fiu, /xie and /xiy), respectively. It is worth 
recalling that the indications of robustness provided by 
the sensitivity coefficients must be taken with caution, 
keeping in mind that this type of analysis is only valid 
locally, i.e., in the neighborhood of the nominal values 
reported in Table 2. 

We next determine the critical points of the system, 
i.e., the points at which system's dynamics undergo abrupt 



changes. We have conducted a bifurcation analysis with 
respect to those parameters that exhibit large sensitivity 
values. Taking Gallp concentration as the output of our 
model, the interval of bistability with respect to a certain 
parameter is delimited by the pair of limit points forming 
the classical S-shaped bifurcation curve. As an example, 
the bifurcation diagram generated by variation of ks is 
shown in Figure 5: the admissible range of variation for the 
parameter ks is ([12.57, 26.62]); outside this interval the 
system loses its bistable behavior. The bifurcation analysis 
can also be performed by allowing simultaneous varia- 
tions of two parameters: in this case, the bistability thresh- 
olds, corresponding to the limit point bifurcations, are 
curves in the parameters plane. For example, in Figure 6, 
where /cy, ks have been chosen as bifurcation parameters, 
we have detected two cusp bifurcation points at (/cy, ks) = 
(4.363,3.835) and (kyjcs) = (166.2,3797.0). Thus, the 
shaded region in Figure 6 identifies a set of parameter 
values within which any value of the pair (/cy, ks) guaran- 
tees bistable behavior (assuming that the other parameters 
take their nominal values). 

Unfortunately, the above methods can efficiently 
evaluate changes in steady state concentrations only for 
low-dimensional parameter variations. In [20], a global 
sensitivity method, named bioSDP, is proposed to evaluate 
the effect of multiple (simultaneous) parameter variations 
on the system's dynamics. More specifically, this method 
can be used to compute some bounds on the maximum 
variation of the equilibrium points induced by changes 
of the parameter values. The computation is based on 
the solution of a dual problem (see Methods for more 
details): given a subset of the state space, say the bioSDP 
algorithm is able to certificate that, for any admissible real- 
ization of the parameter vector p, Q. does not contain any 
equilibrium point. 
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Sensitivity IVIatrix - Time Integral 




Figure 4 Local sensitivity analysis at the steady states of the species involved with respect to all model parameters. The sensitivity matrix 
denotes the normalized sensitivity coefficients of all species in the model across 20 hours. The sensitivity of the parameters k], /t2 (associated in the 
internalization of external galactose and in activation of GalBp protein), ks (basal expression of transcription factor Gal4p), kj, kg, kg (associated in the 
feedback control), /ii3, /ii6 and /ii/ (degradation of GalBp, internalized galactose and Gal2p proteins, respectively), can more influence the 
dynamical behavior of this mechanism more than the other parameters. 



The bioSDP algorithm takes as inputs the admissible 
range of variation of the parameters, defined as a box 
Bp in the parameter space, and an initial outer approxi- 
mation, in the form of a box 5", of the subset Xe of all 
the admissible equilibrium points of system (1) subject to 
p e Bp. 

Then, these outer boundaries are iteratively narrowed 
by applying a bisection algorithm. As a result, the state 



space is partitioned in one or more subsets contain- 
ing all the equilibrium points that fall inside the initial 
search subspace 5". In fact, due to the computational 
burden, the bisection algorithm can only be applied to 
systems of dimension less or equal than three (see for 
example Figure 5 in [21], where a two-dimensional sys- 
tem is analyzed). For higher-order systems like ours, 
the algorithm resorts to a box shrinkage procedure, i.e.. 



25 
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10 12 14 16 18 20 22 24 26 28 30 

ks 

Figure 5 One-parameter bifurcation diagram. The diagram has the classical S-shape in the interval [1 2.57, 26.62], thus the system is bistable for 
values of ks belonging to this interval. (LP, Limit Point). 
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it just tries to reduce the size of the initial box as 
far as possible. 

Due to the above limitations, in our case the method 
proposed in [20] would not be able to distinguish the 
two distinct steady state subsets. To overcome this issue, 
we devise a strategy that leverages the bioSDP algorithm 
but, instead of computing one set containing all the equi- 
librium points, aims to separately compute two distinct 
robust steady state subsets Si and 52, which define the 
boundaries for the variation of Xei and Xe2i respectively. 
Thus, we need two initial outer approximation subsets, let 
us denote them by S'^ and iSj , respectively. 

Guessing two good initial outer approximations would 
in general turn out to be a daunting task for systems of 
nontrivial dimension. In our case, exploiting the previous 
analysis and by virtue of continuity arguments, we sur- 
mise that, for small-enough variations of the parameter 
values, the DA's represent good initial guesses. Thus, we 
let Sf = PiVi, i = 1, 2, where pi > 0 is a scaling factor 
and apply the bioSDP algorithm separately to these two 
initial boxes, with pi = 1; if the algorithm does not find a 
solution, it is re-applied iteratively with different values of 
Pi in a predefined interval [ Pmin> Pmax]. until a solution is 
found. 

Note that, setting as the initial search space for the 
bioSDP algorithm we are focusing the analysis on those 
equilibrium points that belong to a neighborhood of Xet, 
instead of searching for all the equilibrium points. Per- 
forming this analysis separately, first on Xei and then on 
Xe2, enables us to ascertain whether the bistability is pre- 
served against parametric perturbations: the answer is 
affirmative if we are able to compute two disjoint robust 
steady state subsets, i.e.. Si H S2 = 0. If this problem 



is feasible for an assigned parameter box Bp, then we are 
guaranteed that bistability is preserved for all p belonging 
to Bp. 

The initial outer approximations used for our anal- 
ysis are reported in normalized parallel coordinates in 
Figure 2b. Their numerical values are 

S'^ = [ 25.214, 1251.211] x [ 25.214, 1331.199] 

X [ 0.0, 15.991] X [ 0.0, 16.0] x [ 0.175, 12.001] 
x[ 4.515, 96.242] x [0.648, 110.404] 
x[0.0, 12.001] x[0.0, 12.001], 

S^ = [ 52.898, 20852.800] x [ 52.326, 22664.821] 
X [ 4.619, 2701.988] x [ 2.236, 123.599] 
X [ 0.0, 2.80] x[ 0.278, 12.400] x[ 22.269, 6112.385] 
X [ 1.737, 235.601] x [ 1.987, 1283.6] . 

To alleviate the computational burden of the 
procedure, the multi-parametric sensitivity anal- 
ysis has been limited to the parameters subset 
0 := {ki, 1(2, ks, kj, ks, kg, fiu, fiie, fJ.17}, which, according 
to the local sensitivity analysis, have a major influence 
on the system dynamics (see Figure 4). The robustness 
has been evaluated against increasingly larger ranges 
of parameter variations, corresponding to ±2%, ±5%, 
±10%, ±20%, ±30% and ±50%, with respect to the 
nominal values given in Table 2. Figure 2 displays the 
computed robust steady state boxes for the various cases. 
The bistable behavior of the GAL regulatory network is 
guaranteed for parametric variations up to ±20% with 
respect to the nominal parameters value. For such uncer- 
tainty values, indeed, the computed subsets Si and S2 are 
still disjoint, since the intervals of Gsa and G8o,3a are not 
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overlapping. For parametric variations of ±30% or more, 
the intersection of the two subsets is no longer empty (see 
Figures 2, panels g and h); in the latter case, it is no longer 
possible to guarantee that the system preserves bistability 
for all admissible parameter values. 

Discussion and conclusions 

Robustness, intended as the capability to cope with fluc- 
tuations of the molecular species concentrations, caused 
by endogenous and exogenous perturbations, and to pre- 
serve biological functions despite inter-individual vari- 
ability of kinetic parameters, is a key feature of biological 
systems. This fundamental feature poses an important 
challenge when trying to describe biological phenomena 
by means of mechanistic mathematical models: a set of 
differential/algebraic equations which, for some value of 
the parameters, interpolates experimental data, cannot 
be considered a plausible model if it does not possess 
the aforementioned robustness properties. Recognizing 
the power of these arguments as tools for testing novel 
biomodels, and with the aim of supporting the valid- 
ity of our proposed model of the GAL regulatory net- 
work, we have devised an analytical procedure which can 
be exploited to investigate the robustness properties of 
biomodels of bistable biological systems. The procedure 
exploits several complementary methods for the analysis 
of nonlinear quadratic systems (i.e., mass-action mod- 
els), devised both by our group and by other authors, 
and its effectiveness has been demonstrated by applying 
it to thoroughly characterize the robustness of bistability 
for a new model of the galactose metabolism regulatory 
system. 

The procedure consists of three phases: in the first 
phase, the properties of the nominal system (i.e., param- 
eters values are assumed to be certain) are investigated, 
since the first requirement is that the reaction network 
is structurally compatible with the existence of multiple 
equilibrium points. This can be ascertained through the 
use of CRNT, which also allows the computation of a can- 
didate set of parameter values. Subsequently, the second 
stage of the procedure focuses on the analysis of the DA's 
of the equilibrium points, using the method devised in 
[13], since the DA can be regarded as a robustness mea- 
sure against perturbations that push the system away from 
its steady state operative condition. The third phase of 
the procedure consists in the analysis of the robustness of 
the system's bistability with respect to parametric uncer- 
tainty. Traditionally, this analysis is based on sensitivity 
and bifurcations analysis; however, these tools are rather 
limited, due to their inability to take into account multi- 
ple simultaneous parameter variations. To overcome these 
limitations, we have proposed a multi-parametric robust- 
ness analysis strategy: by opportunely leveraging a global 
sensitivity analysis method, and combining it with the 



information provided by our DA's analysis technique, we 
were able to certify the persistence of bistability in the face 
of multiple variations of the uncertain parameters. 

Beyond its specific application for validation of the pro- 
posed model of the GAL regulatory network, the overall 
procedure provides a powerful approach for the analysis 
and validation of any biochemical network model which is 
required to robustly reproduce bistable dynamics, under- 
lying persistent memory, molecular switches and cell dif- 
ferentiation phenomena, without recourse to large-scale 
numerical simulations. 

Methods 

Chemical reaction network theory 

Given a reaction network, the capability of the asso- 
ciated model to exhibit two or more equilibrium 
points depends on the mathematical form of the reac- 
tion rates and on the specific values of the kinetic 
parameters. 

While the characterization of multistability for a generic 
nonlinear system requires an ad hoc mathematical treat- 
ment, in the case of mass-action systems it can be 
performed through a powerful analytical tool, namely 
Chemical Reaction Network Theory (CRNT) [23,24]. 
CRNT links the structure of a biochemical network, 
endowed with mass-action kinetics, to the capability of 
the network to admit multiple positive steady-states. The 
advantage of CRNT is that it provides an immediate 
way to analyze the type of dynamical behavior that one 
can expect from an arbitrarily complex reaction net- 
work, just by inspection of the topology of the associated 
graph. More specifically, CRNT enables us to establish 
the conditions for the existence, multiplicity and stabil- 
ity of fixed points for the associated ODE system, without 
even the need to write down the kinetic equations nor to 
assign values to the kinetic parameters. This point makes 
CRNT especially suitable for dealing with biomolec- 
ular systems, whose parameters are often unknown 
or subject to significant variability among different 
individuals. 

For a complete description of the CRNT, the inter- 
ested reader is referred to the original articles [23,24], 
or to [28] for an introductory overview of the main 
results. Despite the complexity of its theoretical foun- 
dations, the application of the main CRNT results is 
straightforward through the use of the CRNT algorithm, 
which is implemented in the CRNT Toolbox^. Given 
the reaction network that we want to study, the algo- 
rithm establishes whether the associated mass-action 
dynamical system can admit multiple positive steady 
states for some values of the kinetic parameters. In the 
affirmative case, the algorithm also computes a set of 
values of the kinetic parameters for which the system 
is multistable. 
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Analysis of the domains of attraction 

The exact computation of the DAs for a nonlinear system 
is generally a very hard problem to solve, especially for 

systems of medium/high order. The problem of comput- 
ing estimates of the DA's has been studied for a long time 
and several methods, based on Lyapunov functions, were 
originally proposed, e.g. in [29,30]. More recently, Chesi 
has devised novel results concerning DA analysis based 
on Sum of Squares (SOS) representation of polynomials 
and Semi-Definite Programming (SDP), [31,32]. More- 
over, effective examples of the usefulness of SOS/SDP- 
based approaches to elucidate the properties of biological 
systems are provided in [33,34]. 

Topcu et al. in [35,36] have dealt with the topic of esti- 
mating a robust DA in the case of uncertain parameters. 
Compared to the method used in this work, their results 
can deal with a larger class of systems, namely polynomial 
dynamical systems; however, they cast the problem in the 
form of Bilinear Matrix Inequalities (BMI's), whose solu- 
tion is much more demanding than LMI's (used by our 
approach) and, thus, its practical applicability is limited to 
systems of low order with few optimization variables.. 

When dealing with nonlinear quadratic systems, an 
alternative approach to DA analysis, proposed by Amato 
and coworkers in [13,14], can be adopted: this method 
allows one to check whether an assigned box (or, more 
in general, a polytope) in the state space belongs to the 
domain of attraction of a given equilibrium point. Such 
problem can be cast as a LMI feasibility problem [37], 
which can be effectively tackled through effective off-the- 
shelf numerical tools. In what follows we provide a brief 
overview of the main result used in the present work. 

First, recall that a box (or, more generally, a polytope) 
S C can be described as follows 



S =C0nv{x^i),X(2)i ■ ■ ■ >X^r)} 

={xeW^ ■.a/x< l,k=l,2,. 



.,q} 



(4) 



where r and q are suitable integers, X(i) denotes the i-th 
vertex of <S, and conv[ ■] indicates the convex hull of the 
argument. System (2) can be rewritten as 



x(t) = Ax(t) + F(x), 



(5) 



with 



F{x) = 



X^F2 



WfJ 



(6) 



whereof, e R''^''. We can now precisely state the prob- 
lem to be solved. Note that, for the sake of brevity, the 
statement of the problem refers to the zero equilibrium 
point, i.e. the origin of the state space. Nevertheless, it 



is easy to generalize the definition and the entire proce- 
dure to non-zero equilibrium points, via a change of state 
variables, as shown in [13]. 

Problem 1. Assume that each eigenvalue of matrix A in (5) 
has strictly negative real part (i.e., the origin is an asymp- 
totically stable equilibrium point); then, given a box S, 
with the origin of the state space lying in the interior 
of S, establish whether S belongs to the DA of the zero 
equilibrium point. 

The following Theorem provides sufficient conditions 
to solve Problem 1. 

Theorem 1, The box S defined in (4) belongs to the DA 

of the zero equilibrium of system (5) if there exist scalars 
y G (0, 1), c > 0 and a symmetric positive-definite matrix 
Psuch that 



1 yaj,^ 



> 0,k= 1,2,. ..,2« 



yak P/c 
X(i)'^{Plc)X(i) < hi = 1,2,..., 2", 
( X(i)'^Fi \ 

X(i)'^F2 

y{A^P + PA) 



(7a) 
(7b) 

(7c) 



\X(i/F„ / 

+ ^1 ^X(i) f 2 ^^(0 ... F/x(i))P<0, i = l,2,..., 2" 

For a fixed y, conditions (7) constitute a set of LMI's, 
which can be easily solved through off-the-shelf efficient 
numerical tools (e.g., the LMILAB provided in the MAT- 
LAB Robust Control Toolbox [38]). 

In order to find the largest possible estimate of the DA, 
Theorem 1 can be applied iteratively, starting from a small 
initial box Vo, surrounding the equilibrium point, and 
then stretching the box at each iteration along the dif- 
ferent dimensions of the state space, until conditions (7) 
become unfeasible. 

Local sensitivity analysis 

Sensitivity analysis unveils to what extent each parameter 
influences the behavior of a given model and, thus, rep- 
resents a first evaluation of the model's robustness. Our 
sensitivity analysis is conducted according to the method 
illustrated in [27], which is based on the computation of 
the sensitivity coefficients: the sensitivity coefficient s,y is 
defined as the normalized partial derivative of the state 
variable Xi with respect to the parameter pj, that is 



Sij{Xi,Pj, t) 



dXiPj 
dpj Xi ' 



(8) 
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The set of differential equations that constitutes the 
dynamical system is coupled to the equations of the sen- 
sitivity coefficients. This allows computing a numerical 
solution to the whole set of equations, thus simultane- 
ously obtaining both the state variables and the associated 
sensitivity coefficients. 

It is worth recalling that traditional sensitivity analysis 
methods are only valid locally with respect to a particular 
point in the model's parameter space, i.e., in the neigh- 
borhood of a certain parameter set. Another significant 
limitation consists in their capability to consider the sen- 
sitivity of the model with respect to the variations of just 
one parameter at a time; indeed, a model might display 
low sensitivity to single parameter variations, while being 
extremely sensitive to simultaneous multiple parameter 
changes. 

Bifurcations analysis 

Bifurcation analysis is concerned with the study of how 
parameter variations affect the number, type and location 
of attractors, e.g., equilibrium points of a dynamical sys- 
tem. Let us consider a generic nonlinear system, with state 
variables denoted by x, and depending on a parameters 
vector p, 

m =f{xit),p). (9) 

A bifurcation occurs at values of p such that small changes 
of the parameters can dramatically alter the number or 
types of attractors of system (9) [39]. 

The changes in the map of equilibrium points can be 
effectively visualized by using a 'bifurcation diagram', in 
which the steady state value of one state variable is plotted 
as a function of a bifurcation parameter. The calculation of 
bifurcations can be performed through continuation soft- 
wares, like the MatCont package [40], which we have used 
to perform the analysis illustrated in Figures 5 and 6. 

Bifurcation diagrams are powerful tools in order to 
investigate the robustness of nonlinear biomodels in the 
face of parametric uncertainty. However, it is necessary to 
take into account that a) analytical solutions for bifurca- 
tions are only available for low-dimensional models, and 
b) that bifurcation diagrams are practically applicable only 
to study the effect of one or two parameters variation at a 
time. 

Global sensitivity analysis via Infeaslblllty certificates 

In view of the limitations of the approaches presented in 
the previous two sections, while it is possible to employ 
them to obtain preliminary information on the parametric 
robustness of a given model, particular care must be taken 
in drawing any conclusions about global properties of the 
system under investigation. 



To overcome these limitations, we have exploited a 
global sensitivity analysis technique for biochemical net- 
works proposed in [20]. Given the admissible parameters 
variation box, the approach proposed by Waldherr et al. 
allows one to compute an outer approximation, S, of the 
region of the state space that contains all the equilibrium 
points, denoted by Xg. The problem can be formalized as 
follows. 

Problem 2. Given system (2) and a box Bp in the param- 
eter space, compute a box S such that S 5 X^, where 

Xe = {x^W\3^p^Bp:N v(x,p) = 0}. (10) 

Note that, apart from trivial cases, the calculation of 
an analytical form of Xe is practically impossible. More- 
over, computational brute-force approaches are applicable 
only to very low-order systems. Monte-Carlo techniques 
can be applied in the other cases, although they may 
require a large computational effort and guarantee only 
probabilistic results. 

Problem 2 can be effectively solved via the method pro- 
posed by Waldherr et al, which can be formulated in the 
form of the following feasibility problem 

imdx€R",p €W" 
s.t.Nv(x,p) = 0 

(11) 

Pj,mm < Pj < Pj,max, ) = l,...m 
Xi,min ^ Xi < Xi^maxt Z = 1, • • • M, 

where /Jy.min. Pj,max> j = 1, ■ ■ ■ , w, define the admissible 
parameter box Bp, and 

*^i,min» ^/,max; i — Ij • • • j are the 
extremal values of the box S of the state space to be tested 
as a candidate solution to Problem 2. 

The optimization problem (11) is not easy to deal with 
from the computational point of view. However, it can be 
tackled by solving its dual version, that is the problem of 
computing regions of the state space that are guaranteed 
not to contain any steady state for any parameter value 
in Bp. The latter can be relaxed to become a SDP prob- 
lem [41], and solved by means of computationally efficient 
convex optimization tools. For a detailed description of 
this procedure, the reader is referred to [20,21]. In this 
works, the computation of a solution to problem 2 con- 
stitutes the core of an iterative procedure, implemented 
by the bioSDP algorithm: starting from an initial large 
region of the state space, the algorithm tries to compute 
one or more partitions containing X^. The procedure is 
very effective for low-order systems (« < 3), since in this 
case a bisection algorithm can be used for the partitioning. 
For systems of higher order, a box shrinkage procedure 
is employed, which can only return one partition S and, 
therefore, is not useful for analyzing the persistence of 
bistability. 
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In order to solve this problem, we have devised an alter- 
native strategy, which combines the results of the DA's 
analysis with the bioSDP algorithm and has proven to be 
effective in the analysis of our case study. The details of 
this approach have been already reported in the results 
Section 'Local and global analysis confirm robustness of 
the bistable equilibria to parametric uncertainty'. 

Endnotes 

^The CRNT algorithm is implemented in the CRNT 
toolbox, which is freely available at http://www.che.eng. 
ohio-state.edu/~feinberg/crnt/ 
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